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ABSTRACT 

We present a detailed parameter study of collapsing turbulent cloud cores, vary- 
ing the initial density profile and the initial turbulent velocity field. We systematically 
investigate the influence of different initial conditions on the star formation process, 
mainly focusing on the fragmentation, the number of formed stars, and the result- 
ing mass distributions. Our study compares four different density profiles (uniform, 
Bonnor-Ebert type, p oc r -1,5 , and p oc r -2 ), combined with six different supersonic 
turbulent velocity fields (compressive, mixed, and solenoidal, initialised with two dif- 
ferent random seeds each) in three-dimensional simulations using the adaptive-mesh 
refinement, hydrodynamics code FLASH. The simulations show that density profiles 
with fiat cores produce hundreds of low-mass stars, either distributed throughout the 
entire cloud or found in subclusters, depending on the initial turbulence. Concentrated 
density profiles always lead to the formation of one high-mass star in the centre of 
the cloud and, if at all, low-mass stars surrounding the central one. In uniform and 
Bonnor-Ebert type density distributions, compressive initial turbulence leads to local 
collapse about 25% earlier than solenoidal turbulence. However, central collapse in the 
steep power-law profiles is too fast for the turbulence to have any significant influence. 
We conclude that (I) the initial density profile and turbulence mainly determine the 
cloud evolution and the formation of clusters, (II) the initial mass function (IMF) is 
not universal for all setups, and (III) that massive stars are much less likely to form 
in fiat density distributions. The IMFs obtained in the uniform and Bonnor-Ebert 
type density profiles are more consistent with the observed IMF, but shifted to lower 
masses. 

Key words: hydrodynamics - instabilities - stars: formation - stars: massive - 
stars: statistics - turbulence 



1 INTRODUCTION 

The current paradigm of present-day star formation sug- 
gests that stars are born in molecular clouds, permeated by 
supersonic turbulence (lElmegreen & Scalo 2004, Mac Lowl 
|&; Klessen||2004| |Ballesteros-Paredes et al.||2007| ). The cores 



have sizes of a fe w tenths of a p ar sec, are very dense with 



(n) ~ 10 b cm" 3 (Beuther et al. 2007), and m many cases 



they show large line widths, indicating supersonic, turbu- 
lent motions with a power-law spectral velocity distribu- 
tion consistent with P(k) ock ~ 2 (jZuckerman &; Evans|l974[ 
Larsonp981 Heyer & Brunt 2004), and thus steeper than 
the Kolmogorov spectrum of turbulence, P(k) (xk~ 5 ^ 3 . The 
steeper power-law exponent is a result of the compressible 
cascade of interstellar turbulence ( Federrath et al.|2010b ), in 



contrast to the incompressible cascade in Kolmogorov tur- 
bulence. The star-forming regions are observed to be frag- 



mented with a filamentary, fractal-like structure (Scalo 1990 
Men'shchikov et al.|2010| and reference therein) . Very dense 



cores that are supposed to form massive stars have higher 
temperatures (T ~ 20 K) in contrast to less dense clouds 
with 10 K ( |Beuther et al.|20Q7[ [Ward- Thompson et al.|2007 ). 



Despite different fragmentation structures and differ- 
ent local environments, the overall interplay of physical pro- 
cesses that contribute to the formation of stars seems to be 
very robust in producing prestellar cores and finally stars 
with a mass distribution that does not show significant dif- 
ferences in most observed regions of our local Universe. 
This mass distribution can be described by a universal ini- 



tial mass function (IMF) (Scalo 1986 1998 Kroupa 2001 
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Chabrier 2003 ). Only under extreme circumstances, i.e. close 



to the Galactic Centre, may the initial mass function differ 
from the universal one. Whereas |Lockmann et aL| (|2010) find 
that even there star formation is consistent with the canon- 
ical IMF, Bartko et al. (2010) clearly exclude a standard 



IMF in favour of a top-heavy mass function in the Galactic 
Centre stellar disks. 

We know from observations that star formation is a 
complex interplay between a number of physical processes 
and ingredients: gravity, turbulence, rotation, radiation, 
thermodynamics, and magnetic fields. However, to what 
extent the various processes have a dominant impact on 
the evolution in comparison to the initial conditions of the 
molecular cloud is still unclear. Especially the impact of the 
initial conditions on the formation of massive stars, the spa- 
tial distribution of stars, and the mass evolution is unknown. 
Observations reveal that massive stars form early and with 
a tendency to be located at the centre of the cloud, whereas 
stars with lower masses form further out and at later times 
Hillenbrand fc Hartmann||1998| 



( Hillenbrand 



1997 



Fischer 



et al.|1998||de Grijs et"aH20 02 ; Siria nni et al.|2002||Goulier- 
mis et al.|2004||Stolte et al.|2006||Sabbi et al.|2008) . 



Theoretical approaches reproduce consistent star for- 
mation key data with a variety of different numerical meth- 
ods, initial setups, and physical processes (see review by 
Kles sen et al.|[2009| . However, a systematic study of how 
the initial conditions influence the fragmentation process, 
the collapse of the gas into stars, the number of stars, and 
their accretion history is still missing. Especially how the 
formation of massive protostars depends on the interplay be- 
tween initial density profile, turbulence, and accretion model 
needs to be studied systematically. The large variety of exist- 
ing numerical simulations all with different initial conditions 
does not allow for a useful comparison. 
Bate fe Bonnell|(|2005|),|Bate| ( |2009a|b|c ) 



Bate et al. 



(2003), 



Clark et al. 1(2008), 



Bonnell et al.|fl2003[|2004| ), and |Bonnell fe Batej ( |2005| ) used 
uniform density distributions with solenoidal (divergence- 
free), decaying turbulent motions on different cloud scales. 
They use a turbulent power spectrum, P(k) oc /c -2 , consis- 
tent with supersonic turbulence, however, the influence of 
different mixtures of initial modes of the turbulence were 
never investigated. In particular, Bate| ( |2009b ) concluded 
from the similarity of their results with two different initial 
turbulence spectra, P(k)ock~ versus P(/e)oc/c~ 3 , that dif- 
ferent turbulence in general has no major influence on star 
formation. However, both of the investigated spectra in Bate 



(2009b) are steep, such that the turbulence is dominated by 



the few large-scale modes (low k) anyway. Different mix- 
tures of solenoidal and compressive modes of the initial tur- 
bulence are expected to have a much stronger influence on 
star formation, which we show here. Krumholz et al. (2007 
2010) favour concentrated density profiles with p oc r~ ' 



referring to observations of dense cores. Their decaying tur- 
bulent velocities are based on a power spectrum of the form 
P(k) oc k~ 2 , but not specifying the nature of the modes. In 



contrast, Klessen (2001 ) used driven turbulence on different 
scales to create dense cores self-consist ently with a p oc r~ 2 



density profile in the outer region. lOffner et al. (2008) com- 



pared driven and undriven turbulence with an initial flat 
power spectrum for the wave numbers 3 < k ^ 4. Federrath 
et al.| ( |2008[ |2009[ |2010b[ ) investigated purely driven turbu- 
lence with the two limiting mixtures of turbulent modes: 



1) fully solenoidal (divergence- free) and 2) fully compres- 
sive (curl- free), and found significantly different density dis- 
tributions, with three times larger standard deviations of 
the density probability distribution function in the case of 
compressive compared to solenoidal driving (see also the 
follow-up studies by [Schmidt et al.|2009||20Tol|Seifried et ah] 



20TT)|Price et al.||201lf Since such strongly different den- 
sity fields are expected to lead to very different modes of 
star formation, we also investigate here three mixtures of 
the initial turbulence (compressive, mixed, and solenoidal). 
Here, however, we only apply the different turbulent modes 
as an initial condition, not continuously replenishing them 
by driving. 

In this work, we combine four different extreme den- 
sity profiles with different turbulent velocity fields to study 
the influence of the initial conditions on the formation of 
stars. The mass of the cloud is kept constant for all simu- 
lations. We investigate the fragmentation, the time scales, 
and the stellar distributions with a focus on how different 
initial conditions lead to different morphology and statistics 
of prestellar cores and stellar clusters. 

The paper is structured as follows: section [2] describes 
the initial density profiles and the applied turbulent veloc- 
ity fields for the simulations, as well as the numerical key 
parameters, and the usage of sink particles. In addition, a 
theoretical estimate of the accretion rate for the pocr -2 den- 
sity profiles is calculated. In section [3] we present the results 
of the simulations, followed by a discussion in section|4] Here 
we concentrate on the cloud evolution and the global stel- 
lar properties. A detailed investigation of the spatial stellar 
distribution will be published in a separate paper. Finally, 
in section |5] we summarise our results and conclusions. 



2 NUMERICAL METHODS & INITIAL 
CONDITIONS 

2.1 Global Simulation Parameters 

We simulate the collapse of an initially spherical molecu- 
lar cloud with a radius of Ro — 3 x 10 17 cm « 0.097 pc, 
centred in a cubic computational domain of length Lbox = 
8x 10 17 cm. The gas with a mean molecular weight of p — 2.3 
is assumed to be isothermal at a temperature of 20 K. The 
isothermal sound speed is given by 



k B T 
pm p 



0.268 km s 



(1) 



with the Boltzmann constant fe, the temperature T, the 
molecular weight p, and the proton mass m p . For all runs 
the total mass enclosed within this sphere is 100 M©. The 
resulting average density is (p) = 1.76 x 10 -18 g cm -3 or 
(n) = 4.60 x 10 5 cm -3 , leading to a free-fall time 



3tt 



32 G (p) 



(2) 



of 1.58 x 10 12 s or 50.2 kyr. However this global average time 
is not a good measure for the strongly concentrated density 
profiles, where star formation and gravitational collapse oc- 
curs on much shorter time scales. All of the initial spheres 
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Table 1. Physical parameters 
Parameter 



Value 



cloud radius 


Ro 


3 x 10 17 cm « 0.097 pc 


total cloud mass 


Mtot 


100 M 


mean mass density 


(P) 


1.76 x 10- 18 g cm" 3 


mean number density 


(n) 


4.60 x 10 5 cm" 3 


mean molecular weight 


li 


2.3 


temperature 


T 


20 K 


sound speed 


C s 


2.68 x 10 4 cm s" 1 


rms Mach number 


M 


3.28 - 3.64 


mean free-fall time 




5.02 x 10 4 yr 


sound crossing time 


tsc 


7.10 x 10 5 yr 


turbulent crossing time 


Uc 


1.95 - 2.16 x 10 5 yr 


Jeans length 


Aj 


9.26 x 10 3 AU « 0.23 Ro 


Jeans volume 


Vj 


1.39 x 10 51 cm 3 


Jeans mass 


Mj 


1.23 M© 



List of the physical parameters of the runs, which are the same 
for all setups. 

are gravitationally highly unstable. With the Jeans length 



Aj 



G{ P ) 



9264 AU = 0.46 Ro 



(3) 



the Jeans volume, given as a sphere with diameter Aj, 
reads Vj = ttAj/6 and the Jeans mass of this sphere is 
Mj = Vj (p) = 1.23 Mq. The central region inside the Jeans 
volume is called the 'core' in the following. Accounting for 
the different masses inside the Jeans core due to different 
central mass concentrations M(r — Aj/2) — M core , it is 
useful to define the new average density (p core ) and free-fall 
time (tff° re ) for the core region Vj. An overview of all the 
physical parameters is given in table [I] the core values for 
the different density profiles can be seen in table [2] 

The simulations do not include radiative feedback nor 
magnetic fields. The simulated density range justifies an 
isothermal equation of state. However, the missing heating 
effect due to radiation leads to more collapsing regions than 
in non-isothermal simulations. We therefore over-estimate 
the number of formed protostars, and the presented stellar 
statistics should more be understood as a comparison be- 
tween the runs rather than an exact measurement of the 
IMF. 



2.2 Numerical Code 

The simulations were carried out using the astrophysical 



code FLASH (Fryxell et al. 2000) in version 2.5 which 



integrates the hydrodynamic equations with a piecewise- 
parabolic method (PPM) ( |Colella fe Woodward||1984| ). The 
code is parallelised using MPI. The computational domain 
is subdivided into blocks containing a fixed number of cells 
with an adaptive mesh refinement (AMR) technique based 
on the PARAMESH library (Olson et aL|l999| ). 



2.3 Resolution and Sink Particles 

For the main simulations an effective resolution of 4096 3 
cells was used, corresponding to a smallest cell size of Ax ~ 
13 AU. In order to avoid artificial fragmentation, the Jeans 



Table 3. Numerical simulation parameters 
Parameter 



Value 



simulation box size Lb ox 

smallest cell size Ax 
Jeans length resolution 

max. gas density Pmax 

max. number density n max 

sink particle accretion radius r acC r 



5 x 1U X ' cm 
13.06 AU 
^ 6 cells 

g cm" 
6.45 x 10 9 cm" 3 
39.17 AU 



General simulation parameters that are related to the numerical 
resolution. 



length 



Aj 



G Pn 



(4) 



at this effective resolution has to be resolved with at least 
4 grid cells (iTruelove et al.||1997). With sink particles, the 



accretion radius has to be at least 2 grid cells at the highest 
level of refinement in order to fulfil this criterion. In our 
simulations we use an accretion radius of 3 Ax, leading to a 
threshold density p ma x of 



Pmax — 



= 4Gfb = 2 - 46xl °" 14 



As heating of molecular gas begins at 
10 

state seems justified 



g cm . (5) 
density of about 



13 g cm 3 , the assumption of an isothermal equation of 



g., |Lar son 1969). However, we will 



see in the results section that fragmentation is slightly over- 
estimated with the assumption of an isothermal equation of 



state up to these densities (see also, |Krumholz et al.| |2007 
Bate||2009c ) 



We apply the sink particle creation criteria of Federrath 



et al. (2010a) to avoid transient density fluctuations to be 
erroneously turned into sink particles, and thus to avoid 
artificial fragmentation. If the density in a cell on the highest 
level of the adaptive mesh hierarchy exceeds the resolution 
limit, pmax, a spherical control volume with a radius of 3 cells 
at the highest level of refinement (r aC cr ~ 39 AU) around 
that cell is investigated for collapse indicators. An accreting 
Lagrangian sink particle is only formed if the gas in this 
control volume: 

• is converging along all principal axis, x, y, and z, 

• has a central minimum of the gravitational potential, 

• is Jeans- unstable, 

• is gravitationally bound, 

• is not within r aC cr of an existing sink particle. 

The numerical parameters for the sink particles are listed in 
table H 



2.4 Initial Density Profiles 

In the simulations the following four frequently used initial 
density profiles are applied: 

(i) Uniform density profile (Top-hat, TH) 

(ii) Rescaled Bonnor-Ebert sphere (BE) 
(hi) Power- law profile p oc r -1-5 (PL 15) 
(iv) Power- law profile p oc r -2 (PL20). 
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Table 2. Core properties of the different density distributions 



setup M core [Mq] p core [g cm- 3 ] n core [cm- 3 ] t™ re [kyr] t™ Te /t\ 



TH 


1.25 


1.76 x icr 18 


4.60 x 10 5 


49.858 


1.64 


BE 


5.84 


8.33 x 1CT 18 


2.18 x 10 6 


23.061 


2.12 


PL15 


11.12 


1.59 x icr 17 


4.16 x 10 6 


16.707 


2.92 


PL20 


23.02 


3.29 x 10- 17 


8.61 x 10 6 


11.615 


4.20 



Core masses, densities, and free-fall times inside a sphere with diameter of a Jeans length (r core = 
Aj/2 = 7 x 10 16 cm). The free-fall time for the top-hat differs slightly from the theoretical value 
calculated by equation |2j, because the data from this table are the numerical values taken from 
the simulation. 



The profiles are motivated by the following reasonings. 
The TH just reflects the initial conditions in a uniform den- 
sity environment with finite size. Neither initial density per- 
turbations have been established nor does the sphere have a 
developed over-density. The BE profile is motivated by the 
theoretical calculation of an isothermal sphere in hydrostatic 
equilibrium confined by external pressure flEbert|1955[|Bo7r 
nor|19 56). The PL20 profile is the limit of the collapsing BE 
sphere at the end of the evolution process. This density con- 
figuration of a singular isothermal sphere is widely applied 
because its collapse can be described by a self-similar solu- 



10" 



tion with predictable in-fall and evolution properties (Shu 
(1977), section 2.4.4). So far studies with a singular isother- 



mal sphere have only been done without turbulent velocity. 
Finally the PL15 profile, which is an intermediate evolution- 
ary stage of the BE sphere before reaching the PL20 config- 
uration, is motivated by observations. The outer region of 
collapsing clouds is ob served to follow a density distribution 
of the form p oc r" 16 ( |Pirogov||2009| ). 

A comparison of the radial shape for all density profiles 
is shown in figure^ Aj marks the Jeans length for the aver- 
age density (p). These four profiles are extreme setups that 
allow us to follow the influence on the central collapse and 
the fragmentation. 

No initial density fluctuations were applied. The density 
of the surrounding gas in the cubic box around the spheri- 
cal molecular cloud is set to 10 ~ 2 times the gas density at 
the edge of the cloud at r = Ro. The initial temperature 
distribution is a step function with the temperature in the 
cloud envelope 100 times larger than in the inner isothermal 
collapsing cloud, which results in a continuous pressure at 
the boundary r = Rq. 



10" 15 =- 



^ io- ib =- 



<X 10" 17 fcr- 



10- 18 =- 



10" 



-I — I — I I I I 



-I 1 1 — I — I I I I I 

TH 

BE 

PL15 

PL20 



Aj/2 



_i i i i i i 



_L 



1000 



10000 



r [AU] 



Figure 1. Comparison of the four initial density profiles adjusted 
to a total mass of 100 Mq within a radius of 0.1 pc. Aj marks 
the Jeans length for the average density (p). 



f = 6.41 (Ebe rt|1955]|Bonnor|19"56| ). The only free parame- 
ter for this configuration is the central density po . In order to 
better compare this sphere with the other clouds, the cen- 
tral density was first chosen such that the outer radius of 
the sphere yielded the given size of 0.1 pc. Then the density 
at every point was rescaled to fit the total cloud mass of 
Mtot 100 Mq. 



2-4-1 Top-hat 

This density implementation is the simplest profile, describ- 
ing the gas density as a step function 



(p) 

0.01 <p> 



with 



(P) 



Mtc 



V 



for r < Ro, 
for r > Ro 

3M t ot 
AtiRI ' 



(6) 



2.4-2 Rescaled Bonnor-Ebert Sphere 

In hydrostatic equilibrium the critical density profile is de- 
scribed by a Bonnor-Ebert sphere with normalised radius 



2.4-3 Power-law Profiles 

As the power-law profiles p oc r~ p diverge in the centre of 
the cloud, an inner radius has to be denned below which the 
density follows a finite function. In these setups this part of 
the profile is described by a quadratic function: 



ar + c 



B 



(ro) 



for ^ r < n , 
for n ^ r ^ Ro- 



(8) 



The reason for this transition instead of a simple cut-off at 
the inner radius is to avoid artificial numerical effects at 
the boundary n. The value for r± was set to 3 (5) times 
the cell size at the highest level of refinement for p — 1.5 
(p = 2.0). The choice for the values of a and c allow for a 
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continuous transition for the density function value as well 
as for the derivative dp/ dr. For p = 1.5 the two values read 
a = 2.227 x 10" 44 g cm- 5 and c = 1.784 x 10" 14 g cm- 3 , 
the values for p = 2.0 are a = 5.804 x 10 -44 g cm -5 and 
c = 1.107 x 10" 13 g cm -3 . The outer radius Ro was set to 
the radius of the cloud, the constant B scales the density 
profile to a total enclosed mass of M to t = 100 M©. Its value 
depends on the inner radius r±. However, for small radii n, 
which is roughly three orders of magnitude smaller than Ro 
in the numerical setup, B converges to 



lim B 



Mtot(3-p) 1 



47T 



(9) 



Depending on the effective resolution and therefore the pa- 
rameter n, the maximum density changes significantly. 

2.4-4 Power-law Profile p oc r~ 2 and Self- similarity 
Based on the analytic treatment of the collapse of a singular 



isothermal sphere by Shu (1977), the evolution of a density 
profile with the general form 



pM>o) = ^, 



k B T 



(10) 



pm p 

can be described using the dimensionless similarity variable 

r 

x — — , 
c s t 



(ii) 



where G is the gravitational constant. The density distribu- 
tion, the mass accretion rate, and the in-fall velocity can be 
transformed to 



P(r,t) 



a(x) 
4nGt 2 



(12) 

(13) 
(14) 

with ol(x) — x~ 2 dm/dx such that the collapse proceeds in 
a self-similar way. The two basic differential equations that 
have to be solved in order to find the values for a and v read 



MsisM) -^m(x) 
u(r, t) = c s v(x) 



[(x - v) 
[(x-v) 2 -!} 



2 >1e 

1 da 

a dx 



a (x — v) 

x 

a [x — v) 

x 



(x — v) 
(x — v) . 



(15) 



The initial density profile must have the form 

c 2 A _o 



p(r,t = 0) 



4ttG 



(16) 



with A > 2. This equation can be rewritten for the PL20 
density setup as 

p(r, t = 0) = qr~ 2 with q = 5.30 x 10 16 g cm" 1 (17) 

for a total enclosed mass of 100 Mq. The constant A in this 
setup has the value 



A=^«61.9. 



(18) 



Comparing the factor A to the number of Jeans masses in 
the cloud 



Mj 



T 5/2 



6 G 3 /2 p i/2 



(19) 



Nj 



Mj 



it can be rewritten as follows to 



4tt' 



,8/3 r 



N 



2/3 



P 1/3 M t 2 / t 3 



oc N 



2/3 



(20) 



(21) 



62/3 

In order to find the theoretical value for the accretion 



factor mo = m(r = 0,t = 0) equations (15) have to be 
integrated from a large x to a value close to zero. For a 
critical sphere with A — 2 this factor is mo — 0.95, for 
A = 61.9 it reaches a very high value of mo ^421 (see 
figure [2J . This finally gives a theoretical accretion rate of 
^3 



Msis m 



G 



x 10 



" 3 M yr- 1 . 



(22) 



The accretion factor mo can be fitted with a power-law de- 
pendence 

4I.52 



mo oc A 



(23) 



(see right plot in figure [2]) which in turn gives a theoretical 
accretion rate close to a linear dependence on the number 
of Jeans masses 



mo oc Nj 



(24) 



2.5 Initial Turbulence 

2.5.1 Power Spectrum of the Turbulence 

The turbulence is modelled with an initial random velocity 
field, originally created in Fourier space, and transformed 
back into real space. The power spectrum of the modes 
is given by a power-law function in wavenumber space (k 
space) with oc /c -2 , corresponding to Burgers turbulence 
(the value for incompressible, Kolmogorov turbulence would 
be Ek oc k~ 5 ^ 3 in this notation), which is consistent with the 
observed spectrum of interstellar turbulence (e.g., |Larson 
1981||Heyer fe Brunt|2004| . The velocity field is dominated 
by large-scale modes due to the steep power-law exponent, 
—2, with the largest mode having the size of the simulation 
box. Thus, changing the slope of the power spectrum is not 
expected to affect the results significantly (see, Bate 2009b). 
However, the random seed and the mixture of modes of the 
initial turbulence can potentially change the results more 
strongly, which we investigate in this study. Concerning the 
nature of the k modes, compressive (curl-free) are distin- 
guished from solenoidal (divergence-free) ones. The simula- 
tion uses three types of initial fields: pure compressive fields 
(c), pure solenoidal (s), and a natural (random) mixture (m) 
of both. These choices were motivated by the strong differ- 
ences found in driven turbulence simulations using purely 
solenoidal and purely compressive driving of the turbulence 



(Federrath et al.|2008, 2009 2010b). Note however that only 



decaying turbulence with compressive, mixed, and solenoidal 
modes are considered here. For each of these three types, 
two different random velocity seeds are created, leading to 
six different initial velocity fields in total (c-1, c-2, m-1, m- 
2, s-1, s-2), which are combined with the different density 
profiles. 

No overall global rotation is imposed on the cloud. Due 
to the random nature of the turbulence, the net rotation, 
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and the net angular momentum are not strictly zero. The 
ratio of rotational to gravitational energy is of the order of 
a few times 10 -3 . 

2.5.2 Mach numbers 

All setups have supersonic velocities. Due to different den- 
sity concentrations and the resulting different refinement 
structure of the AMR grid, the rms velocities and their Mach 
number 



differ slightly among the different density profiles. Table [4] 
shows the Mach numbers for all the setups which vary from 
M = 3.28 - 3.64 with an average of (M) = 3.44. 



turbulence and density profiles. Table [4] gives an overview 
of the combinations. The BE profiles as well as the PL15 
profiles are combined with all turbulent fields. As the TH 
runs are computationally very expensive, only the turbu- 
lent fields with mixed modes are applied. The PL20 density 
distribution has a very short central free-fall time and is 
expected to collapse and form a massive sink particle be- 
fore the turbulent motions have an important impact on the 
cloud structure. Therefore 3 additional setups with compres- 
sive velocity field c-1 but higher rms Mach numbers (PL20- 
c-lb, PL20-c-lc & PL20-c-ld) were simulated. The velocities 
in PL20-c-lb are twice as high as the ones in PL20-C-1; runs 
PL20-c-lc and PL20-c-ld have velocities 4 and 6 times as 
high as PL20-C-1. The rms Mach numbers are: M c -ih = 6.57, 
Mc-ic = 13.1, Mc-id = 19.7 (see tab.|2}. 



2.5.3 Sound Crossing Time and Turbulence Crossing 
Time 

The sound crossing time through the entire sphere is 

t sc (R ) = 7.10 x 10 5 yr, (26) 

about one to two orders of magnitude higher than the core 
free-fall time for the TH or the PL20 profile, respectively. 
For the supersonic turbulence with an average gas velocity 
of Mach 3.44, the average turbulence crossing time is 

t tc (R ) = 2.06 x 10 5 yr. (27) 

The crossing times for the core region are £g° re = t S c(Aj) = 
1.64 x 10 5 yr and t t c ° re = ttc(Aj) = 4.77 x 10 4 yr, which is 
close to the global free-fall time. 

2.6 Runs 

In order to systematically investigate the influence of the 
initial conditions, we follow a variety of combinations of 



3 RESULTS 

We followed the collapse to a star formation efficiency of 
20%, i.e., until 20% of the initial cloud mass was captured 
in sink particles. The concentrated profile PL20 reached that 
stage quite quickly (~ 11 kyr). The PL15 runs show large 
differences in the simulation time, ranging from 25 — 36 kyr, 
which is similar to the time needed for the BE density setups 
(27 — 35 kyr). The longest time was needed for the TH setup 
with 45 - 48 kyr. Table [5] gives an overview of the total 
simulated time for all setups. Related to the core free-fall 
time, the TH and PL20 profiles just need roughly one t| ore to 
capture 20 Mq in sink particles, whereas the BE runs need 
1.2-1.5 tff re . The longest time was needed by the PL 15 
profiles with 1.4 — 2.1 £^ ore . A comparison of the captured 
mass in sink particles can be seen in figure [3] for all runs. 
The setups with the same density profile are plotted in the 
same line style in order to keep the plot readable. 

During the collapse of the cloud two different gravita- 
tional processes compete with each other. Firstly, the col- 
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Table 4. List of the runs and their main properties 



density 
profile 


turbulent 
modes 


seeds 


name 


effective 
resolution 


M 


total 

^kin 
\E V ot | 


total 

£ therm 
\E V ot | 


core 

Vl^potiyc 


core 

( ^therm 

\ l-Epotl 


TH 


mix 


1 


TH-m-1 


4096 3 


3.3 





075 





047 


0.027 


0.038 


TH 


mix 


2 


TH-m-2 


4096 3 


3.6 





090 





047 


0.111 


0.038 


BE 


compr 


1 


BE-c-1 


4096 3 


3.3 





058 





039 


0.073 


0.028 


BE 


compr 


2 


BE-c-2 


4096 3 


3.6 





073 





039 


0.055 


0.028 


BE 


mix 


1 


BE-m-1 


4096 3 


3.3 





053 





039 


0.018 


0.028 


BE 


mix 


2 


BE-m-2 


4096 3 


3.6 





074 





039 


0.082 


0.028 


BE 


sol 


1 


BE-s-1 


4096 3 


3.3 





055 





039 


0.057 


0.028 


BE 


sol 


2 


BE-s-2 


4096 3 


3.5 





074 





039 


0.072 


0.028 


PL15 


compr 


1 


PL15-C-1 


4096 3 


3.3 





056 





038 


0.067 


0.025 


PL15 


compr 


2 


PL15-C-2 


4096 3 


3.6 





068 





038 


0.042 


0.025 


PL15 


mix 


1 


PL15-m-l 


4096 3 


3.3 





050 





038 


0.013 


0.025 


PL15 


mix 


2 


PL15-m-2 


4096 3 


3.6 





071 





038 


0.072 


0.025 


PL15 


sol 


1 


PL15-S-1 


4096 3 


3.3 





053 





038 


0.052 


0.025 


PL15 


sol 


2 


PL15-S-2 


4096 3 


3.5 





069 





038 


0.061 


0.025 


PL20 


compr 


1 


PL20-C-1 


4096 3 


3.3 





042 





029 


0.046 


0.017 



PL20 compr 
PL20 compr 
PL20 compr 



1 PL20-c-lb 
1 PL20-c-lc 
1 PL20-c-ld 



1024 3 6.6 
1024 3 13.1 
1024 3 19.7 



0.170 0.029 
0.682 0.029 
1.534 0.029 



0.192 0.018 
0.768 0.018 
1.728 0.018 



In order to increase the influence of the turbulence for the PL20 profile three more runs (PL20-C- 
lb, PL20-c-lc, PL20-c-ld) were carried out with the same structure of the velocity field as PL20- 
c-1, but with rescaled absolute values by factors of 2, 4, and 6, leading to rms Mach numbers 
A4 c _ib = 2M c -i, M c -ic = 4A4 c -i, M c .id = 6 M c -i- See table |Al| for resolution details. 



Table 5. Overview of the simulation time and the sink particle 
properties 



Run ^sim tsiva/t ( ^ >Ve ^sim/^ff ^sinks ^max 

[kyr] 



TH-m-1 


48.01 


0.96 


0.96 


311 


0.86 


TH-m-2 


45.46 


0.91 


0.91 


429 


0.74 


BE-c-1 


27.52 


1.19 


0.55 


305 


0.94 


BE-c-2 


27.49 


1.19 


0.55 


331 


0.97 


BE-m-1 


30.05 


1.30 


0.60 


195 


1.42 


BE-m-2 


31.94 


1.39 


0.64 


302 


0.54 


BE-s-1 


30.93 


1.34 


0.62 


234 


1.14 


BE-s-2 


35.86 


1.55 


0.72 


325 


0.51 


PL15-C-1 


25.67 


1.54 


0.51 


194 


8.89 


PL15-C-2 


25.82 


1.55 


0.52 


161 


12.3 


PL15-m-l 


23.77 


1.42 


0.48 


1 


20.0 


PL15-m-2 


31.10 


1.86 


0.62 


308 


6.88 


PL15-S-1 


24.85 


1.49 


0.50 


1 


20.0 


PL15-S-2 


35.96 


2.10 


0.72 


422 


4.50 


PL20-C-1 


10.67 


0.92 


0.21 


1 


20.0 


PL20-c-lb 


10.34 


0.89 


0.21 


2 


20.0 


PL20-c-lc 


9.63 


0.83 


0.19 


12 


17.9 


PL20-c-ld 


11.77 


1.01 


0.24 


34 


13.3 



The time of each simulation is given as the absolute time £ s i m , 
the time in core free-fall times £ s imAff° re 5 and the time in average 
free-fall times £ s imAff- ^sink shows the number of sink particles 
at the end of the run, M max gives the mass of the most massive 
sink particle. 



lapse toward the centre of mass and secondly the collapse of 
dense regions into filaments, induced by the turbulence. The 
different density profiles and turbulent fields lead to different 
cloud evolutions, fragmentation properties, and sink parti- 
cle accretion rates. A column density plot at the end of each 



20 



15 

if io 

5 





5 10 15 20 25 30 35 40 45 50 
t [kyr] 

Figure 3. Comparison of the total mass in sink particles M for 
all simulations. All velocity realisations for one density profile are 
combined in one line style. A detailed discussion of each velocity 
field is given in the analysis section of each of the density profiles. 

simulation is shown in figure [4] and [5] Figure [4] shows the 
column density plots for the density profiles TH, BE, and 
PL15 with the velocity field c-1, c-2, m-1, and m-2, as well as 
PL20-C-1. Each picture row shows simulations with the same 
initial turbulent velocity field, each column belongs to one 
density distribution. In the upper part of figure [5] we show 
the final column density for the BE and PL 15 profile with 
the solenoidal fields. The lower part shows the PL20 pro- 
file with compressive turbulent modes for realisation 1. The 
four different plots belong to different initial kinetic energy 
variations (see table [4|. All simulations show the formation 
of filamentary structures and sink particles. Depending on 
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PL20-C-1 t = 11 kyr 



A^sink = 1 



1 1 1 1 I I I I I 1 1 1 1 I I I I I 1 1 1 1 — I I I I I 

i i i I i i i I i i i I 

lcr 1 io° io 1 io 2 

column density [g cm -2 ] 

Figure 4. Column density plots for the TH, BE, and PL15 setups with velocity profiles c-1, c-2, m-1, and m-2 as well as for PL20-C-1 
at the end of the simulation. The box in all cases spans 0.13 pc in both x and y direction. Each picture row corresponds to one velocity 
field, each column to a density profile. All setups show filamentary structures but differently spread in the box. Only the TH density 
runs form distinct subclusters. 



the initial density profile, the turbulent field, and the re- 
sulting total simulation time, the position of the filaments 
as well as the number of sink particles and their spatial 
distribution vary significantly. The TH profiles in figure ^] 
show locally disconnected filaments and subclusters of sink 
particles. The BE profiles also form many sink particles in 
extended filaments, but much more centrally concentrated 
and in stronger connected filaments. The initial mass con- 
centration and the resulting faster central collapse suppress 



the formation of completely disconnected subclusters. The 
PL 15 density profile shows in many cases a similar cloud 
evolution as the BE setups. However, the total number of 
sink particles varies strongly with different velocity realisa- 
tions and the sink particles are located closer to the centre of 
mass. The influence of different initial kinetic energies of the 
turbulent motions can be seen in the PL20 setups. Higher 
velocities lead to much stronger substructures within the 
same simulated time. 



Importance of the Initial Conditions for Star Formation I: Cloud Evolution and Morphology 9 




10" 



10° 



10 1 



10 2 



column density [g cm 



Figure 5. Column density plots for the BE and PL15 setups with velocity profiles s-1 and s-2 (upper part) as well as for the PL20 setup 
with turbulent field c-1, c-lb, c-lc, and c-ld (lower part). The box in all cases spans 0.13 pc in both x and y direction. 



A time evolution for turbulent field m-2 and the den- 
sity profiles TH, BE, and PL15 is shown in figure [6] Each 
row shows the column density at the same simulation time. 
The columns correspond to the different density profiles. 
The much slower central collapse in the TH case allows 
the formation of two distinct over-dense regions, shown at 
t — 22 kyr. At that time the BE profile has formed a few 
stars along the long main filament. The PL 15 profile has al- 
ready formed more than 50 sink particles very close to each 
other that interact very strongly and disturb the central fil- 
amentary structure. 3 kyr later the BE sphere formed more 
stars mainly along the outer arms of the main filament. Al- 
though the number of sink particles is larger than in the 
PL 15 case at the previous time snapshot and the total mass 
captured in sink particles is roughly comparable, the cluster 
is not dominated by the gravitational attraction and Af-body 
dynamics of the stars. The initial gas structure remains un- 
perturbed. Another 3 kyr later the TH profile eventually 
developed collapsing regions in completely disconnected ar- 



eas. By that time the BE cluster begins to show dynami- 
cal interactions. In the last time snapshot the overall cloud 
structure as well as SFE and the number of sink particles is 
comparable for the BE and the PL15 case. 

Concerning the formation of sink particles, a clear dis- 
tinction between the power-law profiles and the profiles with 
a flat core has to be made. The power-law profiles with their 
high density core form a sink particle very early due to the 
fast collapse of the central region. In the PL20 profile and in 
two of the PL15 profiles, this particle remains the only par- 
ticle formed in the entire simulation time. PL15 runs with 
more than one sink particle form them with a large time gap 
after filamentary structures have formed and collapse. In the 
BE and TH profiles this central particle does not exist, and 
all particles form in filaments. This different behaviour can 
be seen in the mass evolution (figure[3|. The runs with PL15 
profile form a sink in the centre right after the start. The 
mass therefore evolves similarly at the beginning. For the 
BE sphere and the uniform density distribution, the differ- 



10 Girichidis et al 



t = 22 kyr I BE-m-2 



t = 22 kyr I PL15-m-2 
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t 25 kyr 
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Figure 6. Column density plots for the TH, BE, and PL15 profile with the velocity field m-2 in a box of 0.13 pc in x and y direction. 
The TH-m-2 clearly develops two subclusters by the end of the simulation. The BE and PL15 runs show a similar general cloud structure 
that is dominated by central collapse. In the BE case the flatter initial density forms sink particles far away from the centre, whereas in 
the PL15 run the cluster is more compact. 
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ent realisations of the turbulence lead to different filamen- 
tary structures and thus influence the point in time when 
sink particles are created. Therefore, the mass evolution of 
the different simulations show large offsets (figure [3|. 

In general, all setups result in high total accretion rates 
onto the sink particles of M ~ 1 - 2 x 1CT 3 M yr -1 . 
Only PL15-m-2 and PL15-S-2 (see detailed discussion be- 
low) show somewhat smaller values of the accretion rate. 
The fluctuations around the mean value strongly depend on 
the number of particles, their positions, and the resulting 
particle-particle interactions as well as accretion shielding 
effects. The PL20 as well as two PL 15 runs only form one 
sink particle and show a very smooth accretion rate with 
small fluctuations. The accretion rates in the TH and BE 
profiles are influenced by the particle movements but as the 
clusters are not that compact the interactions are less in- 
tense. The overall similarity of the accretion rates can also 
be seen in the similar slope of the mass function in the upper 
panel of figure [3] 

3.1 Analysis of the TH Profile 

The uniform density distribution has much less mass within 
the core region compared to the concentrated profiles (see 
table [5]), and its core free-fall time is longer. The initial 
supersonic velocity field has time to develop significant 
over-densities before the global collapse becomes dominant. 
Therefore the evolution of the cloud at the beginning of the 
simulation is dominated by the turbulent motion rather than 
the central collapse. The turbulence crossing time and the 
free-fall time of the core are similar (tt? re /t| ore =1.64) which 
leads to the formation of over-dense regions all over the sim- 
ulation box. These over-dense regions are very massive and 
evolve to locally collapsing filaments in which the first sink 
particles form. Filaments that are close to each other merge 
into sub-cores in which subclusters build up, before the cen- 
tral collapse sets in. After roughly one free-fall time, 20% of 
the mass is collapsed into sink particles. 

The accretion rate for every single sink particle is a 
strongly varying function with time. However, the global 
rate for the sum of all sink particles quickly reaches a satu- 
rated value of M ~ 10 -3 M© yr _1 (figure m), which can also 
be seen in the comparable slope of the total sink particle 
mass as a function of time. The number of sink particles is 
noticeably higher for TH-m-2. 

The mass distribution of the sink particles follows an 



20 



overall shape similar to the universal IMF (e.g. Kroupa 2001 
Chabrier 2003), but shifted to lower masses by a factor of 
about 10 (see figure [8|. A comparison with analytic models 



of the IMF (e g., |Padoan & Nordlund||2002l jHennebelle ~& 
Chabrier|2008 ) is planned in a future contribution. Here the 
main conclusion is that the formation of massive stars is 
very unlikely in a cloud with 100 Mq and a uniform density 
distribution. 

Since refinement is initiated in a very space-filling fash- 
ion for the uniform density distribution of the TH runs and 
thus computational cost became prohibitive, we only ran 
mixed turbulence runs with two different seeds. It should 
be noted, however, that the influence of the different mix- 
tures (compressive versus solenoidal) of the turbulence has 
the biggest influence on the evolution and structure of the 
forming clusters and subclusters in the TH profiles, because 
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Figure 7. Sink particle evolution for the TH runs. Once sink 
particles have started to form the total accretion rate (lower plot) 
quickly reaches a value around M ~ 10 -3 Mq yr -1 , fluctuating 
by a factor of roughly 2. Therefore the evolution of the mass 
captured in sink particles as a function of time looks very similar 
for both runs (upper plot), just shifted by 3 — 4 kyr. 



TH profiles provide the most time for the turbulence to in- 
fluence the cloud structure before the global collapse sets in. 
This will be addressed in a separate paper. 



3.2 Analysis of the BE Profile 

Here the cloud evolution at the beginning is similar to the 
collapse of the TH core. The turbulence can form strong fila- 
ments spread over large regions of the domain. However, the 
different radial mass distribution leads to low-mass filaments 
in the outer regions. This results in a stronger central accel- 
eration, which causes the filaments to merge near the centre 
of mass. The formation of large subclusters is suppressed 
compared to the case of the uniform density distribution. 
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Figure 8. The mass distribution of the sink particles for the TH 
setup has an overall shape similar to the uniform IMF ( |Kroupa| 
|2001||Chabrier|2003| ), but shifted to lower masses. 



By far, most of the sink particles, which are roughly as nu- 
merous as in the TH simulations, are formed in the core 
region. The time evolution of the cloud for different turbu- 
lent modes with the same random velocities can be seen in 
figure [9] The compressive modes lead to sink particle for- 
mation about 25% earlier than the mixed and solenoidal 
modes. 

The time evolution of the global sink particle properties 
are shown in figure [l0| Although the random seed strongly 
determines the location and orientation of the filaments, the 
particle formation between BE-c-1 and BE-c-2 is almost in- 
distinguishable. In the case of mixed and solenoidal modes 
the choice of the random seed significantly changes the time 
at which sink particles form. However, after the creation of 
sink particles has set in, the particle production rate with 
time as well as the total mass accretion rate is quite similar 
for all runs, not reflecting the structure of the initial turbu- 
lence at all. Only the BE-s-2 setup needs some more time un- 
til it reaches the asymptotic value of M ~ 2 x 10 -3 Mq y _1 . 
However, the accretion rate of individual sink particles varies 
strongly with time. 



The mass distribution of the sink particles (figure 1 1 ) 
also shows a typical IMF structure like the TH runs, also 
shifted to much lower masses. This leads to the conclusion 
that the stronger central density concentration and the re- 
sulting stronger in-fall properties are still way too inefficient 
in forming massive stars. 

3.3 Analysis of the PL15 Profile 

From the very beginning of the simulation, the PL 15 pro- 
files show a considerably different evolution compared to the 
TH and BE profiles. Due to the strong mass concentration, 
the first sink particle forms close to the centre very early, 
after roughly 1 kyr ^0.06 t^ ore . The formation of this sink 
particle is not influenced by extended filaments, because the 



20 



15 



10 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 m; 1 1 1 1 1 i,i i'i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

BE-c-1 / 

BE-c-2 ./ / 

E" BE-m-l / //, : 

BE-m-2 / /// '' 

BE-s-1 

BE-s-2 





300 
250 
200 
150 
100 
50 


0.01 



(3 0.001 



0.0001 



1 I I I I I I i \\ \ I I I I I i if I i Ti I I I I 1 1 1 I I I I I I I 1 1 1 I I I I I I I 1 1 1 I I I I I I IT 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ! 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



u 



ITm?i i i i il'i in I I I i 



I I I I I I I 1 1 I I I I I I I I 1 1 I I I I I I I I 1 1 I I I I I I I I 1 1 I I I I I I I I 1 1 I I I I I I LL| 



; 7 



ill 



10 



15 



20 



25 
t [kyr 



30 



35 



40 



Figure 10. Sink particle evolution in the BE runs. The upper 
plot shows the total mass captured in sink particles. The com- 
pressive fields form sink particles first, the mixed and solenoidal 
velocity fields a few kyr later. After the formation of the first 
sink particle the accretion rate (lower plot) approaches a value of 
~ 2 x 10 — 3 Mq y — 1 , independent of the initial turbulent field. 
The number of sink particles also shows a similar evolution for 
all setups (central plot). 



formation time of filaments is much larger than the time for 
central collapse. The central particle has a high and smooth 
accretion rate in all PL 15 runs, which allows it to grow to the 
most massive sink particle in the simulation, while filaments 



in the outer regions start to form later (figure 12). Whether 
secondary sink particles form strongly depends on the ran- 
dom seed of the turbulence, as well as on the nature of the 
modes. All simulations with compressive modes lead to the 
formation of many sink particles in the filaments. On the 
other hand, mixed and solenoidal modes lead to either one 
(PL15-m-l & PL15-S-1) or a few hundred particles (PL15- 
m-2 & PL15-S-2). A possible explanation for this dichotomy 
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Figure 9. BE column density plots for the BE density profile and three different turbulent fields. The columns show snapshots with c-1, 
m-1, and s-1 velocities (from left to right) for the same physical time. The box shown spans 0.13 pc in x and y direction. 
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Figure 11. IMF for the BE setups. For all turbulent setups 
the IMF looks very similar. The distribution function mainly fol- 
lows the general shape of a uniform IMF, but with lower average 
masses. 



could be the influence of tidal forces, which can suppress 
the growth of the initial perturbations induced by the turbu- 
lence. In a density profile steeper than r _1 (see appendix |b|, 
tidal forces start to shear radial density fluctuations apart, 
thus reducing the chance of initial perturbations to grow by 
self-gravity. For the BE profile the central region of the cloud 
has a shallower density profile than r _1 , the PL15 profile a 
slightly steeper one. Over-dense regions that can marginally 
grow in the BE profile may be sheared apart in the corre- 
sponding PL 15 profile with the same velocity field. However, 
the turbulence is supersonic and the density power-law ex- 
ponent is not far from the critical one. This is why different 
locations and strengths of converging and diverging regions 
of the velocity field may easily overcome the shearing effect 
and cause the big differences between PL15-m-l/PL15-s-l 
and PL15-m-2/PL15-s-2. Indeed, an analysis of the density- 
weighted divergence of the initial velocity fields shows that 
seed 2 produces stronger compressions in regions of high 
density than seed 1. Taken together with the fact that frag- 
mentation into multiple objects always occurs for the purely 
compressive fields, this shows the importance of compressive 
modes for triggering the formation of dense fragments. 

In the first 10 kyr, the evolution of all PL15 simulations 
is quite similar. During that time all simulations have only 
formed one central sink particle. As soon as other sink par- 
ticles form, the situation changes significantly. In the case 
of multiple sink particles, their particle-particle interactions 
in the stellar cluster disturb the central in-fall and redirect 
the central gas velocities. 

Although the total number of sink particles as a func- 
tion of time is similar for PL15-C-1, PL15-C-2, PL15-m- 
2, and PL15-S-2, their spatial distribution differs between 
the runs with compressive velocity field (PL15-C-1, PL 15- 
c-2) and the runs PL15-m-2 and PL15-S-2 with mixed and 
solenoidal fields. In the former, the sink particles are located 



in filaments much farther away from the centre, resulting in 
weaker particle-particle interactions and allowing the parti- 
cles to remain located in their dense parental filament. The 
runs PL15-m-2 and PL15-S-2 are dominated by the in-fall of 
less centrally located and hence less massive filaments. The 
local gravitational collapse inside these filaments is there- 
fore delayed until the filament approaches the dense core. 
Sink particles show much lower mean separations which in- 
creases the strength and impact of particle-particle interac- 
tions. The induced cluster dynamics reduces the total mass 
accretion rate because individual sinks stop accreting if they 
are kicked out of the dense gas regions. This effect can also 



be seen in the IMF (figure 13). PL15-m-2 and PL15-S-2 have 



many more sink particles, but the final mass of the central 
one is lower than in the runs with compressive fields (see ta- 
ble]^]). Hence, the accretion onto the central object is starved 
by the fragmentation around it ( Peters et al.||2010a| ) 



3.4 Analysis of the PL20 Profile 

For the PL20 density profile with the compressive turbulent 
field, only one sink particle was created already after 0.13 kyr 
which is only 0.012 t^ ore . As this velocity field is the most 
likely one to form more than one sink particle, the other 
turbulence realisations are not simulated entirely. This den- 
sity profile is gravitationally too unstable for the turbulence 
to have an impact on the density evolution and the frag- 
mentation of the gas sphere within a core free-fall time. As 
the turbulence crossing time is about 20 times longer than 
the core free-fall time, the small influence of the turbulence 
is expected. The accretion rates for all realisations of this 



setup are very similar (figure 14). Therefore only the setup 
with compressive mode 1 (PL20-C-1) was simulated up to a 
star formation efficiency of 20%. In conclusion, a p(r) oc r~ 2 
density profile does not reproduce a realistic IMF but helps 
to form massive stars. 

In order to investigate the threshold turbulent energy 
that is needed to cause other regions to fragment and col- 
lapse besides the central region, three additional PL20 pro- 
files with higher velocities were investigated (see tab.[4j). The 
turbulence in PL20-c-lb, with twice as high velocities than 
our standard PL20 run, is still not strong enough to signif- 
icantly alter the result. There is still only one sink particle 
created, accreting mass at a very high rate. For PL20-c-lc 
with velocities four times as high as in PL20-C-1 (A4 = 13.1), 
the situation changes. The stronger turbulence leads to the 
formation of other sink particles apart from the central one. 
However, the central particle in this run still contains 90% 
of the mass (M = 18 M©) at the end of the simulation, and 
the second most massive particle is more than one order of 
magnitude less massive. Similar results are obtained from 
PL20-c-ld with a Mach number of M = 19.7. More sink 
particles form, but still the central star is the most massive 
one with M = 13 M©. 



4 DISCUSSION 

Our results clearly show that diverse initial conditions lead 
to completely different cloud structures and collapse sce- 
narios. However, the strong dependence of the simulation 
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Figure 12. PL15 particle evolution. The upper plot shows the 
total mass captured in sink particles. Apart from the m-2 and s-2 
velocity field, the mass evolution is very similar. This can also 
be seen in the lower plot, showing the accretion rate. In case of 
more sink particles, the accretion rate varies much more strongly 
with time. This is due to strong particle-particle interactions in 
the compact stellar cluster. If the cloud fragments and collapses in 
different regions the number of stars follows similar curves (central 
plot). However, the total number of particles differs much more 
than in other density setups. 



outcome on the initial conditions may be moderated by dif- 
ferent input physics like radiation or magnetic fields and the 
effects due to rotation. 

Our simulations indicate that massive stars can form 
without the aid of radiation and magnetic fields just from 
choosing centrally concentrated density profiles. In contrast, 
our isothermal cloud setups with flat density distributions 
fail to produce massive stars. We note that this result could 
change significantly if more massive clouds with more Jeans 
masses are used. We find in our simulations of isothermal 
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Figure 13. IMF for the PL15 runs. All runs form one very mas- 
sive sink particle, which is by far the most massive one in the 
cluster, indicated by the single peak around 10 Mq. The continu- 
ous set of low-mass particles below the mass gap shows similarities 
with the universal IMF, again shifted to almost 10 times lower 
masses as in the TH and BE profiles. 



gas that clouds with an initially uniform density distribu- 
tion tend to overproduce low-mass proto-stars and have dif- 
ficulty forming sufficient numbers of high-mass objects. This 
is in qualitative agreement with the simulations done by 
andlBate h Bonnelll (120051). They used a 



Bate et al. 



(2003) 



uniform density distribution and solenoidal velocity fields, 
which seems to represent the conditions that inevitably lead 
to a large number of low-mass objects. This is also consistent 



with the calculations by Offner et al. (2008), Klessen et al 



(1998 ), Kless en fe Burkert| fl2000p |Klessen et al.| ( |2000a| >7 
|Klessen ~ p)01| ), and |Heitsch et aL] ( |2001| ), who tested the 
influence of driven and decaying turbulence in a uniform 
density box. In order to suppress fragmentation and/or en- 
hance the formation of massive stars in flat density profiles, 
more physics may help, which is addressed in three differ- 
ent approaches, namely radiation feedback, magnetic fields 
and stellar collisions. Concerning the first process, |Kratter| 



& Matzner (2006) derived an analytical model to address 
the fragmentation process in massive discs. Indeed, Bate 



(2009bp, |Krumholz et al.| ( |2009| ), |Peters et al] ( |2Q10a|c|b| ) 

found reduced fragmentation in simulations. However, ra- 
diative feedback does not suppress fragmentation entirely. 
Alternatively, magnetic fields tend to reduce fragmentation. 
Hennebelle & Teyssier| ( [20081 ), |Ziegler| ([2005]), and |Biirzle| 



et al. (2010) investigated the influence of magnetic fields in 



low-mass cores, Banerjee et al. (2009), Peters et al. (2011), 



and Hennebelle et al. (2011) noted reduced fragmentation 



in high mass cores. But again, the fragmentation is not fully 
suppressed. The formation of massive stars by stellar colli- 



sions was proposed by 


Zinnecker & Yorke (2007). However, 


Baumgardt h Klessen 


(2010) showed that under realistic 



cloud conditions the contribution of stellar collisions can be 
neglected. 
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Figure 14. Comparison of sink particle accretion rate for the 
PL20 profile with different turbulent velocity fields. Note that the 
accretion rate is plotted in linear scale in order to see the small 
differences between the runs with different turbulence. In all cases 
only one sink particle is created over the whole simulation time. 



In earlier studies, the discussion about cloud fragmen- 
tation and the formation of massive stars is strongly focused 
on the physical processes, not taking heed of the importance 
of initial conditions. As the time scale for star formation is 



of the order of a few dynamical times ( Ballesteros-Paredes 


et al.|1999; Elmegreen 2000 


|Hartmann et al.|2001||Mac Low 


& Klessen|2004| Elmegreen 


2007 ) , the star- forming core has 



only little time to interact with the surrounding medium 
The boundary and initial conditions are therefore decisive 
key properties for the collapse scenario and the star forma- 
tion outcome. To fully understand the formation of a star 
cluster therefore requires knowledge of both the initial con- 
ditions for the cluster- forming cloud core (density profile, 
temperature, turbulent velocity content) as well as the time- 
dependent boundary conditions (as the core is connected to 
the overall turbulent cloud environment and may grow in 
mass by accumulation of gas at the stagnation points of 
larger-scale convergent flows). 

Given the sensitivity of the dynamical evolution on the 
choice of the initial density profile, it is of pivotal impor- 
tance to seek guidance from observations. On small scales 
(<C 1 pc) the observed cores clearly deviate from a uniform 
density (e.g., Pirogov 2009; Ko nyves et al.pOlO Bontemps 
et al.|2010| ). The outer regions of molecular cloud cores can 
be described by a power-law with p oc r -1 ' 5 . In the centre of 
a dense core, however, the approach of a power-law function 
seems to be inconsistent with observations, which identify 
the central region of the core to be flat flMotte et aL]|1998| 



Ward-Thomp son et al.||1999| ). Starless cores may often be 
fitted with a critical Bonnor-Ebert sphere, cores with stars 



are often better fitted with super-critical ones (Teixeira et al. 
Kandori et al.|2005||Kirk et al.|2005) . |Krumholz et al 



2005 



(20071 12010) use very similar setups to our PL 15 density pro- 
file, emphasising the importance of radiative feedback for the 
formation of massive stars. In this density profile the cen- 



tral region inevitably determines the collapse time scale and 
the formation of the first proto-stellar object. Our current 
analysis indicates that following a power-law profile to very 
small radii (< 10 3 AU) introduces a bias towards forming a 
massive central object without much fragmentation around 
it. Adding radiative feedback does not change the outcome 
significantly in view of the very short central collapse time 
scales. 

We can also look at the way ISM turbulence is treated 



in other numerical studies. 


Bate et al. (2003), Bate & Bon- 


nell (2005), Bate 


(2009a|b|c 


), Bonnell et al. ( 2003 2004 ), and 


Bonnell & Bate ( 


2005), for example, always used divergence- 



free, de caying turbulent fields. |Clark, Glover, &; Kl essen 



( [2008| , |Clark, Bonnell, fe Klessenj fl2008| ), |Offner et~ 



(2008), and Krumholz et al. (2007) do not specify the nature 
of the modes they select for their turbulence. As our results 
show that compressive, decaying modes lead to significantly 
earlier collapse and more elongated, shocked structures in 
the flat density profiles (TH and BE) than purely solenoidal 
turbulence, this is an important aspect of the star forma- 
tion process that deserves further consideration. A system- 
atic study of different modes of the turbulence was done 



by|Federrath et al.| (|2008| |2009| |2010b| ) and |Seifried et al. 
fl2011| ), but' in a periodic box with driven turbulence and 
without gravity. These studies find the expected trend that 
compressive modes initiate faster collapse and higher accre- 
tion rates than purely solenoidal turbulence. However, the 
influence of the different modes is stronger in driven turbu- 
lence with self-gravity than in the decaying turbulence runs 
analysed here. Since dense cores are typically embedded in 
large-scale, turbulent molecular clouds, an effective driving 
of the internal turbulence from outside the core is expected 
Klessen et al.|2000b| |Federrath et al.||2010b ). 



(e.g., 



5 SUMMARY AND CONCLUSIONS 

We performed a parameter study of the fragmentation prop- 
erties of collapsing isothermal gas cores with different initial 
conditions. We combined four different density profiles (uni- 
form, Bonnor-Ebert type, p oc r -1-5 , and p oc r -2 ) and 
six different turbulent, decaying velocity fields (compres- 
sive, mixed, and solenoidal, each with two different random 
seeds). For these simulations we neglected radiation, mag- 
netic fields, and initial rotation, in order to study the direct 
influence of the initial density profile and the character of 
the turbulence. The cloud evolution as well as the star for- 
mation and their properties were examined. Here we list our 
main conclusions: 

The density profile strongly determines the number of 
formed stars, the onset of star formation, the stellar mass 
distribution (IMF), and the spatial stellar distribution. 

• Flat profiles (uniform density and Bonnor-Ebert pro- 
files) produce many sink particles in elongated filaments. 
The formation of sink particles starts after slightly more 
than half of a core free-fall time for the uniform cloud and 
after roughly one core free-fall time for the Bonnor-Ebert 
setups. The runs with initially uniform density produce sub- 
clusters in merging filaments in outer regions of the cloud. 
Even the relatively weak mass concentration in the centre 
of the Bonnor-Ebert setups suppresses the formation of sub- 
clusters. Both density profiles show an initial mass function 
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with the high-mass end consistent with the Salpeter slope. 
In the case of initial compressive velocity fields, star for- 
mation sets in 25% earlier than in the solenoidal case. The 
mixed turbulent fields are in between the two extreme cases. 

• The p oc r -1-5 profiles always form one sink particle in 
the centre of the cloud at an early stage. This sink particle 
accretes gas at rate of ~ 10 -3 M© y _1 and grows to the most 
massive particle by far. The formation of unstable filaments 
depends sensitively on the initial turbulent field. The for- 
mation of additional sink particles only occurs after a time 
delay of ~ 0.3 ts. The mass distribution of these sink par- 
ticles shows a high-mass slope consistent with the Salpeter 
slope, but has a wide gap between this mass continuum and 
the central massive star of almost an order of magnitude in 
mass. The spatial distribution shows a compact structure 
around the centre of the cloud and no subclustering. The 
column density of the filamentary structure looks extremely 
similar for a p oc r -1 ' 5 run and the corresponding Bonnor- 
Ebert run with the same turbulent field, not reflecting the 
significantly different stellar properties. 

• The p oc r~ 2 density profile quickly leads to the for- 
mation of one single, central sink particle. The formation 
of other stars is strongly inhibited due to the rapid col- 
lapse compared to the time scale for filament formation. In 
this scenario further star formation can only be triggered by 
higher Mach numbers of the turbulence, if the ratio of tur- 
bulent energy to gravitational energy is increased to about 
unity. 

The realisation of the turbulent velocity field has a ma- 
jor impact in the different morphology of the filamentary 
structure, their orientation, and shape. 

• In the uniform density profile the random seed of the ve- 
locity determines the position of filaments from which stars 
form, and thus the location of the stellar subclusters. In ad- 
dition, the number of sink particles generally depends on the 
random seed of the turbulence. Similar results are obtained 
for the BE profile. 

• The p oc r -1 ' 5 profile, which marks the transition be- 
tween one central massive sink particle and many low mass 
ones, is very sensitive to the random seed. Different realisa- 
tions may switch between one single star and several hun- 
dred. The formation time and location of the central, first 
sink particle, however, is not influenced by the random seed. 

• The p oc r~ 2 setups are not noticeably influenced by the 
turbulence. The short collapse time of the core compared to 
the turbulence crossing time does not allow for turbulence 
to strongly influence the evolution. 

Our results suggest that massive stars predominantly 
form out of highly unstable cloud cores which are either 
strongly centrally concentrated or much more massive than 
modelled here, allowing stars to accrete form a larger mass 
reservoir. The density configuration with p oc r -1-5 seems to 
be the most sensitive one concerning the number of collaps- 
ing fragments for different turbulent velocities. 

Overall we conclude that the choice of the initial density 
profile is an extremely important, perhaps even the most im- 
portant parameter determining the fragmentation behaviour 
of high-mass proto-stellar cores. Choosing an ideal simpli- 
fied density profile strongly preordains the subsequent star 
cluster properties. This implies that the effects of different 



physical processes can only be reliably compared if the ini- 
tial density profile is the same. In realistic star formation 
simulations, the formation of these cores needs to be taken 
into account and cores need to be formed self-consistently 
from larger clouds. 
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APPENDIX A: RESOLUTION STUDY 

We test the influence of the numerical effective resolution of 
the code on the collapse by simulating the different cloud 
setups with different resolutions. The order of the numeri- 
cal resolution was chosen such that the isothermal approx- 
imation for the equation of state is appropriate and the 
computational effort is acceptable. The different resolutions 
have acronyms corresponding to the maximum refinement 
level (RL) in the code: Z max = 7 (RL07), Z max = 8 (RL08), 
/ max = 10 (RL10) and / max = 12 (RL12). Due to the dif- 
ferent sizes of the smallest cell, the maximum gas density 
before creating sink particles, as well as the accretion ra- 
dius vary. A comparison of the parameters can be seen in 
table [AH 

As the computational time for the BE and the TH pro- 
file are very large (i.e., more than an order of magnitude 
larger than for the PL20 profile, because of the quite space- 
filling refinement in the evolution of these profiles), these 
setups have only been compared in an early evolutionary 
stage. The highly concentrated PL20 cloud has been inves- 
tigated in more detail: for a longer evolution time, for more 
different resolutions and analytically. 



Al BE Profile 

Due to the flat inner core of this profile, refinement is initi- 
ated in a rather large volume of the core, which makes the 
computational effort for this profile much larger than for the 
other profiles, and thus the resolution test was done only for 
a short simulation time. In figure |A1| we compare the to- 
tal accretion rate, M, and the number of sink particles N 
of the Bonnor-Ebert profiles BE-c-1 and BE-s-1 for resolu- 
tions RL10 and RL12. The accretion rates are comparable 
and give roughly the same star formation efficiency with 
time. However, the number of particles varies significantly 
with resolution. This is expected, since we use an isother- 
mal equation of state, which does not introduce a physical 
length scale or density threshold to the problem, i.e., the 
problem remains scale- free. Changes in the equation of state, 
in particular if the gas becomes optically thick, will break 
the scale- free collapse (e.g., Jappsen et al.|[2005 Krumholz 
eTaLp007| [Bate||2009ct . 



A2 PL20 Profile 

For the concentrated density profile with p oc r~ 2 and the 
turbulence profile c-1, detailed simulations were run for four 
different maximum refinement levels: RL07, RL08, RL10, 
RL12. In all cases only one sink particle was created in the 
centre of the cloud after a few steps of hydrodynamical evo- 
lution. 

The results for the PL20 runs can be seen in figure |A2| 
The accretion rate onto the protostar M does not differ sig- 
nificantly, resulting in the same slope of the mass M as a 
function of time. The different evolution of the accretion 
rate at the very beginning of the simulation is due to the 



different geometrical setup conditions (see sec. 2.4.3). The 
larger size of the smallest cell for lower refinement levels re- 
sults in a much coarser density distribution in the centre of 
the cloud and needs more evolution time in order to develop 
a sink particle with constant accretion rate. The theoretical 
value for the accretion rate fits the simulated values very well 
(see sec. 2.4.4). The comparison with a simulation without 



turbulent velocities only shows minor differences. 



APPENDIX B: TIDAL FORCES 

The tidal acceleration in a spherically symmetric setup at 
distance r from the centre with an enclosed mass M is given 
by 

M 



atidai(r) = G 



r±Ar) 2 '' 



(Bl) 



where G is the gravitational constant and Ar <C r. The en- 
closed mass can then considered to be constant within the 
variation Ar. Given a density profile of the form p(r) oc r~ p 



yields a mass function M(r) 
tion scales as 



r 3-p 



and the tidal accelera- 



atidai(r) oc r 1_p . 
The derivative with respect to r, 

dtttidal , 



dr 



-(r) oc (1 — p) r 



(B2) 



(B3) 



changes sign at p = 1. For p < 1, atidai increases with radius 
(datidai/dr > 0) and therefore compresses material at radius 
r. For p > 1, datidai/dr < and shears condensations apart. 
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Table Al. Main simulation parameters for different resolutions 



refinement 


eff. res. 


Ax [AU] 


facer [AU] 


Pmax [g Cm 


" 3 1 


Wmax [cm 3 ] 


RL07 


512 3 


104.4 


313.3 


3.85 x 10- 


16 


1.01 x 10 8 


RL08 


1024 3 


52.2 


156.7 


1.54 x 10- 


15 


4.03 x 10 9 


RL10 


4096 3 


13.1 


39.2 


2.46 x 10- 


14 


6.45 x 10 9 


RL12 


16384 3 


3.3 


9.8 


3.94 x 10- 


13 


1.03 x 10 11 



Main simulation parameters for different effective resolutions. The accretion radius of the sink 
particles r aC cr is set to 3 times the minimum cell size Ax. 
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Al. Comparison of the Bonnor-Ebert profiles BE-c-1 and BE-s-1 for resolutions RL10 and RL12. The mass accretion as a function 
(left plot) shows only small differences between the two resolutions. The number of sink particles, however, differs strongly (right 
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Figure A2. Resolution comparison for the PL20 runs with turbulence field c-1. After an initial evolution time the accretion rates approach 
the same value for all setups. The differences at the beginning of the simulation are due to different maximum central resolutions. The 
flattened density function at the centre of the box is much shallower for lower resolutions resulting in larger times for a central collapse. 



